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Abstract 

The effect of final state interactions on the elliptic flow coefficient measured in relativistic heavy 
ion collisions is investigated within the DWEF formalism established by Miller and Cramer [Phys. 
Rev. Lett. 94, 102302 (2005), J. Phys. G34, 703 (2007), Phys. Rev. C78, 054905 (2008)]. It is 
found that the optical potential previously found to give the best fit of particle multiplicity and 
Hanbury Brown and Twiss radii in RHIC events has a moderate effect on the resulting elliptic 
flow coefficient v%. This indicates that final state interactions should be taken into account to 
confidently predict V2 to better than ~ 20% accuracy. 
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I. INTRODUCTION 



Hydro dynamic models have been quite successful at describing single particle hadron 
spectra measured in experiments at the Relativistic Heavy Ion Collider (RHIC) 
Recently, progress has been made in going beyond simple ideal hydrodynamic models and 
using viscous hydrodynamics to determine the value of various transport coefficients in the 
RHIC "fireball," e.g. the shear viscosity . Of particular interest is the azimuthal 



anisotropy of produced particles. The existence o 
best indications of thermalization at RHIcJiol. \ll 
connected to the viscosity of the medium |6j. 



: strong elliptic flow at RHIC is one of the 
H . [l3 |. and its precise value is sensitively 



Not every aspect of these hydrodynamic models, however, is completely physically jus- 
tified. Due to uncertainty in especially early and late time dynamics in the evolution of a 
heavy ion collision, the initial and final conditions for hydrodynamic evolution are treated 
somewhat simplistically. In particular, the so-called Cooper- Frye freezeout algorithm [l4] is 
typically implemented to describe how the medium transitions from hydrodynamic behavior 
to the essentially free particles that enter the detectors. Using this procedure, the medium, 
once it reaches some defined freezeout temperature, instantaneously "freezes out" from a 
hydrodynamic fluid to completely free particles that stream to the detectors. Knowing that 
this is not entirely physical, it makes sense to investigate how much this treatment might 
affect the conclusions reached from using these hydrodynamic models. 

To that end, Miller, Cramer et al. investigated the effect of introducing some final state 
interactions to this freezeout behavior in the form of a one-body optical potential that the 



15 
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The 



emitted particles interact with as they "fight" to escape the medium 
main motivation behind this previous work was the notorious inability of hydrodynamic 
models to fit two-particle correlation data (so-called Hanbury Brown and Twiss radii) while 

n n 

simultaneously fitting single-particle data [18(. (For a review see Ref. [19[.) Thus only 
multiplicity and HBT radii data were calculated, while the present work extends that to 
now investigate the effect of these final state interactions on measured elliptic flow. For 
simplicity only pions — the dominant hadron produced in a RHIC event — are considered. 

The distorted-wave emission-function (DWEF) formalism is briefly reviewed in Sec. [Til 
with the calculation of v-i in this formalism described in Sec. IHII Sections IIVI and IVl contain 
the results and conclusions. Details of the calculation are included in Appendix |A] for 
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those interested, along with a semi-analytic derivation of a simple test case used to test the 
numerics, in Appendix [B] 



II. DWEF FORMALISM 



The DWEF formalism was established in Ref . [15j and described extensively in Ref . 16] , 
and the relevant parts are briefly summarized as follows. 

The main quantity that we are interested in is the detected particle momentum spectra 

E*£- = -£L-= td±xS{p,x) (1) 



d 3 p dY d 2 p 

from which vi is defined as the second moment in the azimuthal momentum angle 

fd4 P cos(20 p )^. 

v 2 = (cos(20 p )) = f ^ dN (2) 

J a<pp dY d 2 p - 

with p the momentum of the detected particle and Y the particle rapidity. <p p is the angle 
of the particle momentum with respect to the collision plane. 

S(p, x) is known as the emission function. In conventional hydrodynamical models, it 
resides in space on a freezeout hypersurface defined by a surface of constant temperature (or 
other thermodynamic quantity) in the hydrodynamic simulation, with a momentum distri- 
bution at each point on the surface given by the appropriate equilibrium (or off-equilibrium 
in the viscous case) distribution at that given temperature. This freezeout hypersurface 
represents the surface of last scattering, from which free particles are emitted and travel 
directly to the detector. 

Here, instead of running a full hydrodynamic simulation, we follow Ref. 20| and use 



an analytical parametrization of the freezeout surface, similar to one typically found in 
numerical hydro simulations, but with tunable parameters. In addition it is allowed to be 
a more general volume, with finite width in all dimensions, rather than the infinitely thin 
surface obtained in a conventional Cooper-Frye prescription. 

Secondly, instead of freely streaming, particles that are emitted from this surface are then 
made to interact with an optical potential representing interactions with the medium from 
which the particles are escaping. 



Explicitly we have (see Ref. 16j for details) 



cosh ri -n\ 1 -(T-n } ) 2 M± p(b) . , 

S(p,x) = - / — r^e^ 7 - e 2Ar* — W~\x)\ 2 . 3 

' (2tt) 3 V27T At 2 e ( P -u-^)/T _ i^p v j\ v ) 



p is the asymptotic pion momentum, and M± = a/p± + m \- As usual, instead of Cartesian 
coordinates (t, x,y, z), we use the set (r,x,y,r]) or (r, 6,0, 77) with 

i] = arctanh(z/t) 



Vt 2 



b = ^x 2 + y 2 (4) 
= arctan(?//:r). 
b= (6,0). 

z is the beam direction, with the xz (0 = 0, tt) plane the reaction plane. 

ipp{x) are the aforementioned distorted waves (as opposed to plane waves appropriate 
in the absence of interactions). They obey the equation of motion 

d 2 

di 2 

for pions interacting with optical potential 



V 2 -^-U(h)-ml)^-\x)=0 (5) 



U(b) = -(w + w 2 p 2 ) p(b). (6) 

Note that although the medium is time-dependent in principle, for simplicity the optical 
potential is taken here to be time-independent and can be interpreted as a time-averaged 
quantity. 

Whereas in the original DWEF formalism only a rotationally symmetric transverse den- 
sity p(b) was needed (corresponding to central collisions), here we are interested in azimuthal 
anisotropy and so we need to consider a more general form. Specifically we take the modified 



Woods- Saxon profile from Ref. [20|] 



(exp[(-l)^l + 1)2 
p(b) = \ y[K ' a - J , (7) 

(exp[(6,/^ + ^-l)^] + l) 2 



-'■ u ' s 



with R ws = \/^(R 2 + R 2 ). Thus lines of constant density in the transverse plane form 



ellipses with semimajor to semiminor axis ratio 



Lastly we must specify the fluid velocity u, for which we again defer to Ref. [20|. It is 
parametrized using a transverse fluid rapidity r] t (b) 

^{x) = (cosh r] cosh rj t , sinh^cos0b, sinh?^ sin 0;,, sinhr/cosh^). (8) 
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The transverse direction is taken to be perpendicular to lines of constant density. It can be 



20] 



shown that the angle of such a fluid velocity, (fib, obeys 

R 2 

<fi b (<fi) =tan- 1 (-|tan0). (9) 
K y 

The transverse fluid rapidity r] t (b) is first taken to have the same elliptic symmetry as 

QcosW) 2 j_ (6 shi,/,)) 2 rj^^ 



the density, increasing linearly with the "radial" coordinate b = \ j 
added to this is a term proportional to cos(20) representing the amount of elliptic flow built 
up before freezeout 



r}t{b) = r) f b> 



'cos z 



Rl 



+ 



sin 



R 2 



1 + a 2 cos(20)). 



(10) 



The momentum in these coordinates takes the form 



p^ = (Mj_ cosh Y, p± cos (fi p , p±_ sin (fi p , M± sinh Y) . 



We choose to focus on data at midrapidity, Y = 0, and so 



p ■ u = Mj_ cosh 7] cosh rj t — pj_ sinh rj t cos( 



(11) 



(12) 



In all, then, the parameters involved in this model are: Ai], At, r , fi n , T, w Q , w 2 , R x , R y , 
a ws , r]f, and a 2 . We are interested in the effect of an optical potential like the one found to 
give the best fit in Ref. [lfj and so we will keep all of these parameters fixed to those best-fit 
values, and only adjust ^ and a 2 to give reasonable results for non-central collisions. 

It should be noted that the formalism developed is not strictly correct when the optical 
potential is comp 



fit values of Ref. 



ex. (See the discussion in Ref. [171].) We therefore also investigate the best 
ItT ] for a vanishing imaginary part of the optical potential. 



III. CALCULATING v 2 

This section outlines how the calculations are carried out. A set of coupled differential 
equations must be solved numerically to obtain the wavefunctions ifip~\ and then a five- 
dimensional integral must be performed (two of which can be done analytically with suitable 
approximations. ) 
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A. The Wavefunctions ifjp (x) 

ipp satisfies Eq. (jSJ). Since U(b) is independent of t and z, we can write 

^ p ~\x) = e^e^^ib), (13) 

and Eq. (jSj) becomes 

(Vl-U(b)+pl)4-\b) = 0, (14) 

or 

(w + l I b ^w- u{h)+pl ) 4 ~ >{h) = - (15) 

Decomposing ^ ■* and U(b) into angular moments 

oo 

^ _) (b)= £ / m (p,6)(-2) m e im ^), (16) 

m=— oo 

tf(b)=£tf„(&)e** (17) 

n 

results in 

e im0 e -im^ = q ^g) 

So the term in brackets vanishes identically for each m, and we must solve a set of coupled 
differential equations. In practice, every f m above a certain m max is set to zero, and a finite 
set of coupled equations is solved numerically. 

The boundary conditions are the same as for the cylindrically symmetric case — far outside 
the medium one should have a canonically normalized plane wave plus an outgoing wave, 
i.e. 

f m (b > R ws ) = J m {p b) + T m H£\p b) (19) 

with J m and Bessel functions and Hankel functions of the first kind, respectively. 

Details of this calculation can be found in Appendix |A] The program used to calculate 
the wavefunctions was tested in part by comparing to a semi-analytic solution described in 
Appendix [B] 



E 



d 2 1 d _ m 2 2 \ X^n f 
dtf + bdb~l r+P± ) fm ~^ nfn 



i n e im 
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B. Integration 



Once the wavefunctions are found, a five-dimensional integral must be performed: 

/ c% cos(2 P ) / d 4 x S(p, x) 
V2 fd<j) p fd*xS( P ,x) ' 1 } 

The t integral can be done analytically 

-(t-Tq) 2 

rdre = V27rr Ar. (21) 



The 77 integral can also be done analytically with the following approximations (as in 



Ref. 



id) 



— r\ 1 cosh 77 

e 2A i 2 w t^ 1 e Ar >' 2 (22) 



-, Jmax 

~ P (~P-u+^)j/T ( 9 o\ 

e (p-u-^)/T _ I C ' 

i=i 

where the Bose-Einstein distribution is approximated by a sum over Boltzmann distributions 
truncated at some j max , and so 

d V coshr? e" cos W^+^co ShJ?t ) = 2 ^ + ^ cosh ^^J (24) 

Finally, then, for the numerator we have 
/ #s cos(2«/^5fex) 



2 r Mj_ 1 ^ 
— — — > 1 e t 
2vr 3 ^ 



x I rf 2 6 p(b)/ ra (p, 6) £(p, 6) e^"^ + cosh ^ 

x / c% cos (20 p ) e - i ( m - n )^eTP^ sinh ('") cos (^-^), (25) 



and similarly for the denominator. The final three integrals are done numerically 
More details of this part of the calculation can also be found in Appendix |A] 

IV. RESULTS 

As previously mentioned, we would like to determine the effect of adding final state 
interactions to hydrodynamic fits. To gain insight into this, we consider an emission function 
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TABLE I: Best fit parameter sets. The top line (Fit 1) is a general fit 16|] while the bottom line 
(Fit 2) is from a fit where Im(w 2 ) is held at 0.0001 [l7]. 





T 


rj f 


At 


Rws aws 


w 


w 2 


T() 


Ar) 






(MeV) 




(fm/c) 


(fm) (fin) 


(fm- 2 ) 




(fm/c) 




(MeV) 


Fit 1: 


156.58 


1.310 


2.0731 


11.867 1.277 


0.0693 


0.856+i0.116 


9.04 


1.047 


139.57 


Fit 2: 


121 


1.05 





11.7 1.11 


0.495 


0.762+i0.0001 


9.20 


70.7 


139.57 



with parameter values taken from Refs. 16|, ll7| , which give the best description of the single 
particle data in general, and also with the imaginary part of the optical potential held at 
zero (see Table [H Also note that in both fits the chemical potential was fixed at the pion 
mass) . 

We must make alterations to this central collision model to approximate a more peripheral 
collision. The results for a central collision do not unambiguously imply what a peripheral 
collision will look like without appealing to a particular model for the dynamics of the 
system. We therefore choose reasonable parameters to approximately represent a collision 
with impact parameter ~ 7 fm, and then see how the resulting v 2 depends on the strength 
of the optical potential. In principle one could vary all the parameters and do a separate fit 
of all the relevant experimental data (multiplicity, HBT radii, v 2 , etc.) for each of various 
collision centralities. However, the computing time to do so would be prohibitive, and here 
we are most interested in investigating only the effect of the interactions, so we proceed as 
follows. 

First, as in Ref. 



161 ]. we scale down R ws , a ws , and t by the number of participants to 
the 1/3 power, with N part taken from the Glauber model (with the same parameters used 
in Ref. [6J) for an impact parameter of and 7 fm (N part = 377.5 and 171.544). Specifically 



0.7688-R ws . Then we adjust the ratio tt- such that the spatial eccentricity 



(y 2 ) - (x 2 ) 



(26) 



(y 2 ) + (x 2 ) 

has a value of 0.035. This is a reasonable value corresponding to the spatial eccentricity at 
freezeout of hydrodynamic fits of peripheral collisions with impact parameter ~ 7 fm. Note 
that the brackets in Eq. [26] indicate a spatial average with weight given by Eq. [7J while the 
spatial eccentricity in hydrodynamic simulations are typically given with respect to , e.g., 
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(a) 



(b) 
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FIG. 1: (Color online) Calculated v% as a function of momentum with 02 = (a) and 02 = 0.11, 0.10 
(Fit 1, 2) (b). Points with error bars are experimental data for pions at 20-30% centrality from 
the STAR Collaboration [2l|. 



energy density. We nevertheless keep the eccentricity from Eq. [26] fixed at this value with 
an understanding that it is only a rough but still realistic guide to the shape. 

Lastly we must specify how much elliptic fluid flow is built up in earlier stages of the 
collision, represented by the value of (recall Eq. [TOl . First we set 02 = and see what v% 
is generated by interactions with the optical potential in the absence of significant elliptic 
fluid flow (Fig. QJa)). The calculated elliptic flow coefficient t>2 is plotted as a function of 
momentum, along with the relevant experimental data. (Note that p in our calculation is 
the momentum of an asymptotically free pion detected far outside the medium, not the 
momentum of a particle as it is emitted inside the medium, and can therefore be compared 
directly to experiment.) Although we are only able to calculate up to a limited momentum, 
it is clear that final state interactions alone do not generate an appreciable value for u 2 for 
either the general best-fit parameters (Fit 1) or those with a vanishing imaginary part of 
the optical potential (Fit 2). 

Next we increase 02 such that the experimental value for v% is roughly obtained (Fig.[T](b)). 
A value of 02 = 0.11 was required for the parameters from Fit 1, while 02 = 0.10 was sufficient 
to bring the emission function from Fit 2 into the physical regime. One can see that the 
optical potential has a small but non-negligible effect — it decreases v 2 on the order of 10-25% 
of its zero-interaction value with a slightly smaller effect as momentum increases. 
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V. CONCLUSION 



Final state interactions in the DWEF model were found to have a small, though not 
entirely insignificant effect on the elliptic flow coefficient t> 2 . This is in addition to the indirect 
effect of adding final state interactions. For example, adding an optical potential changes 
other observables such as the multiplicity, which would alter parameters in a hydrodynamic 
fit such as freezeout temperature, which would then in turn have an effect on the calculated 
value of t>2. 

The precise size of these effects in general can only be determined with a better under- 
standing of the model fits (e.g. Fit 1 versus Fit 2) in addition to a more detailed analysis — a 
full parameter search using all the relevant experimental data, or perhaps even by adding 
final state interactions directly into hydrodynamic simulations (i.e. a hydrodynamic after- 



burner in the vein of, e.g., Refs 



23 



24 



251]). It is reasonable, however, to conclude 



that final state interactions can affect the calculated value of t; 2 by as much as ~ 20% 
(in agreement with other investigations of final state interactions, e.g., Ref. [l|), and so 
must be properly taken into account to have confidence in the quantitative predictions of 
hydrodynamic simulations at that level of precision. 
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APPENDIX A: DETAILS OF THE NUMERICAL IMPLEMENTATION 

A program was written in C++, making use of the GNU Scientific Library (GSL) version 
1.9, to do the calculation of v%, as detailed here. 

The integral over the azimuthal angle of the pion momentum, p is done as a sum using a 
simple trapezoid rule. This is because for each different value of 4> p , a new set of differential 
equations must be solved. This also allows for the numerator and denominator of Eq. [2D] to 
be solved simultaneously, with just a factor of cos(20 p ) multiplied to the numerator when 
adding terms to the sum. 
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For each term in the sum, then, first the wavefunctions ipp ^ are obtained. They obey a 
set of coupled differential equations of the form 

(w + 11 - £ + pl ) fm - £ u » ^ = ( A1 ) 

for all integers m. This set is truncated, since large m moments (f m for m > pj_R ws ) 
contribute little to the wavefunction. Therefore, all f m for m greater than some m max are 
set to zero, leaving a finite (2m max + 1) number of coupled ordinary differential equations. 
These are solved by calling a GSL solver. Using an embedded Runge-Kutta-Fehlberg method 
seemed to give the best performance. For these solutions, Eq. [7] is integrated numerically to 
find the moments U n . This is done with the GSL adaptive integration routine for oscillatory 
functions. 

To match to the proper boundary conditions, one must find (2m max + 1) linearly inde- 
pendent solutions to this set of equations and take the correct linear combination of these 
solutions that matches the desired boundary conditions. The straightforward choice for 
these linearly independent solutions is to sequentially solve for the case where only one of 
the partial waves is non-zero near the origin. For example, for the n'th solution let: 

fmiP bmin — ) ^m,n 

P 

f m {b m in) = -^&m,n (A2) 

and then solve the set of differential equations up to some arbitrarily large b max far outside 
the potential. We can then match each partial wave in this n th solution to the form: 

J m (pb) + B m>n H£\pb). (A3) 

The final wavefunction is then given by the linear combination of these solutions that 
matches the form of Eq. [L_H] at b max : 

f m (b) = ^2C n f m>n {b). (A4) 

■n 

This part of the program was tested with the trivial case of zero optical potential, in 
addition to comparing to a separately written program that calculates only the cylindri- 
cally symmetric well as to the results of the semi-analytical test case described in 

Appendix [Bl 
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Once these wavefunctions are obtained and stored in memory, the integral over b and 
in Eq. [25] can be performed in addition to the sum over Boltzmann factors. The integrations 
are done with two GSL adaptive integration routines, one embedded in the other. The sum 
is done inside the argument of the integrals. 

APPENDIX B: SEMI- ANALYTIC TEST CASE 

To test the numerics, the case of a pion moving through an elliptically-shaped step- 
function potential was solved (semi-) analytically making use of elliptic coordinates. This 
can be compared to the case of a ws — > (see Sec. HT|) . 

We want to solve Eq. [5] with U(h) an elliptically shaped step function — a finite potential 
inside an ellipse in the transverse plane, with zero potential outside. 

It is useful to change to elliptic (cylindrical) coordinates, denoted u and v. Think of u as 
a 'radial' coordinate that runs from to oo and v as an 'angular' coordinate that runs from 



Note the major and minor axes of the resulting confocal ellipses are reversed from the shape 
of the density used in the main calculation (which is larger in the y direction). This is to 
maintain consistency with the conventional definition of elliptic coordinates. At the end one 
can simply take p — > (<f) p + it) to match the usual convention in RHIC papers. 
Consider the case 



to 2vr 



x = a cosh(u) cos(v) 




(Bl) 




(B2) 



The sharp boundary at u = u is an ellipse with major and minor axes 



R x = a cosh(wo) 



R y = a sinh(ii ). 



(B3) 



In this coordinate system the Laplacian is 




(B4) 
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and so Eq. (jSJ) becomes 



(sinh 2 (u) + sm 2 (v)) \du 2 dv 



d 2 d 2 
+ 



U(u) +p 2 



V P (b) = o 



(B5) 



or equivalently 



with 



d 2 d 2 
— + 2q{u) cosh(2u) + — - 2q(u) cos(2t;) 
ou z ov z 



q(u) = j(p 2 ~ U(u)) . 



i>Jb) = 



(B6) 



(B7) 



On the inside of the potential and on the outside separately, q{u) does not depend on u and 
these cases can be solved with separation of variables and the solutions patched together at 
u = Uq. Let 



Qout 



4 



(p 2 ~ Uo) 



-P 2 . 



(B8) 



Start by expanding if) p (b) in terms of so-called elliptic sines and cosines of the 'angular 



variable v. They are solutions of 'Mathieu's equation' 

d 2 



dv" 1 



+ 2q cos(2w) C(a, q,v) = a C(a, q, v). 



(B9) 



The general solutions are called 'Mathieu functions,' usually denoted C(a, q, v) for solutions 
even in the coordinate v and S(a, q, v) for odd. Demanding periodicity of the variable v 
allows only certain discreet eigenvalues a (denoted here a n for the even functions and (3 n 
for the odd functions). This (complete) set of periodic solutions is commonly called elliptic 
sines and elliptic cosines: 



C(a n ,q,v) = ce n (v,q) 
S(/3 n ,q,v) = se n (v,q). 



(BIO) 



The general solution of Eq. (1B6P can be written in terms of these elliptic sines and cosines: 

oo 

^p( b ) = Yl \fcn{u)ce n (v,q) + f Sn (u)se n (v,q)] . (Bll) 



n=0 
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Plugging this in to Eq. (1B6j) gives 

d 2 



du 2 
du 2 



a,. 



+ 2q cosh(2w 
+ 2q cos!i(2m) - f3 n 



fsM 



-- 
0. 



(B12) 
(B13) 



This is called the modified Mathieu equation, which can be obtained from Eq. (IB9I) by 
replacing v — > (i u). Note that the eigenvalues are different for the functions corresponding 
to ce n and se n (f Cn and f Sn above, respectively). The general solution is then the same as for 
the original Mathieu equation, analytically continued with v — > (i u), though typically they 
are organized by boundary conditions analogous to Bessel and Neumann functions (denoted 
Je n (u,q), Ne n (u, q), etc.) Q: 



f Cn (u) = C Cn Je n (u, q) + S Cn Ne n (u, q) (B14) 
fs n ( u ) = C Sn Jo n (u } q) + S Sn No n (u, q). (B15) 

Note that there are many different sets of so-called Mathieu functions, each being a 
complete orthogonal basis. Replacing with q out results in a different basis, and there are 
separate sets of modified Mathieu functions corresponding to the eigenvalues of the elliptic 
sines and elliptic cosines (a n and f3 n ). 

By requiring continuity at the u = line segment one finds that the general solution 
inside the potential is: 

ip™{u, v) = [Ce™ Je n (u, q in )ce n (v, q in ) + Co*" Jo n (u, q in )se n (v, q in )] (B16) 
with undetermined coefficients Ce^ n , Co™. 



28] 



Outside, we write the solution as the sum of a plane wave and an outgoing wave 

YA -Je n (u, qout) + Ce° n ut He£\u, q out ) 
n ^ Pn ' 

+ (^-Jo n (u, q ou t) + Co^Ho^iu, qout^j se n (v, q ou t)se n ((p p , q ou t)}, (B17) 

where the H's are analogous to Hankel functions 

He£\u, q) = Je n {u, q) + i Ne n {u, q) (B18) 
HoW(u,q) = Jo n (u,q)+iNo n {u,q) (B19) 
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1 1 ^ 



Pn 7T J Q 



and the plane wave coefficients p n and s n are 

dv e ip - x ce n (v,q out ) (B20) 

1 1 f 2w 

- = - dve^ x se n (v,q out ). (B21) 

Sn K Jo 

The coefficients Ce° n ut and Co° u *, along with the analogous 'inside' coefficients are determined 
by matching boundary conditions. 

To match at the u = u boundary project the 'inside' angular functions 
(e.g. ce n (v,q in )) in terms of the 'outside' ones (e.g. ce n (v,q out )). 

oo 

cej(v, q in ) = ^2 Bj n ce n (v, q out ) (B22) 

n=0 

oo 

sej(v, q in ) = ^2 B jn se n(v, q ut), (B23) 



n=0 

with 

-27T 



1 f ZW 

B in = - dvcej(v,q in ) ce n (v,q out ) (B24) 
n Jo 

1 f 2n 

B in = ~ dvsej(v,q in ) se n (v,q out ). (B25) 
7T J 

Then the 'inside' wave functions are 

% n = \° e f Je i( u > qin "> B jn ce n(v, qout) + Cof Joj(u, q in ) B s jn se n (v, q out )] . (B26) 

The coefficients (Ce™, Co^ n , Ce° n ut , Co° n ut ) can then be determined by demanding that ip 
and its gradient be continuous at u = u , which gives the following relations: 

^2 CefJej(u , qin) B j n = —Je n (u , q ut)ce n ((j) p , q out ) + Ce° n ut He (1 \uo, q ut)ce n ((p p , q out ) 

(B27) 

^2 CofJojiuo, qin) B ] n = — Jo n (u , q ou t) se n ((p p , q ou t) + Co^Ho^Xuo, qout)se n ((frp, q ou t) 

(B28) 

^2 CefJe'jiuo, qin) B j n = —Je' n (u , q ou t)ce n ((j) p , q out ) + Ce™*#e /(1) (u , q ou t)ce n ((j) p , q out ) 

(B29) 

^2 CofJo'^Uo, q in ) B j n = —Jo' n (uo, q ut)se n ((f) p , q ou t) + Co° n ut Ho {1) (u , q ou t)se n ((f) p , q Q 



{out) 



(B30) 
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The plane wave coefficients (p n , s n ) as well as the coefficients from the projection 
(Bj n , Bj n ) must be solved numerically. In addition, to compare to the f m in the main 
calculation, the resulting wavefunctions are integrated to project out the usual angular mo- 
ments. Hence the description as a "semi- analytical" test case. In fact, this implementation 
(done in Mathematica) saves no time over the original numerical version, but it does provide 
an independent check. 
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